################### PROJECT- JCR REPLICATION - figures_and_table_main_02 ##########################################
# This R file plots Table 2, Figure 7 and 8  in the main article
# Those figures and table serves the analysis purpose
# Below are brief introduction for those figures and table 
# Table 2 - The Unconditional Impact of Natural Disasters on Two Conflict Outcomes
# Figure 7 - The conditional impact of natural disasters on two conflict outcomes, by territorial setting: 90% and 95% confidence intervals
# Figure 8 - The conditional effect of natural disasters on two conflict outcomes: a three-way interaction between natural disasters, territorial setting, and rebel Control
# Last updated: 16th May 2025, by Wangyin Zhao
##############################################################################################################

##############################################################################################################
# Outline of Script
# A. Preparation
#    - Set working directory
#    - Load required packages
# B. Import Data
# C. Table 2
# D. Figure 7
# E. Figure 8
##############################################################################################################

#### A. Preparation ####
# Set working directory to Replication_Files folder
setwd("./Replication_Files")

# Import relevant packages
packages <- c("tidyverse", "sf", "tidyverse", "multcomp", "fixest", "gridExtra")
lapply(packages, library, character.only = TRUE)
########

#### B. Import data ####
# To generate those figures and table, the first step is to run the data processing files and to generate the analysis data 
# Upload dataset1_oldschool.csv
dataset1 <-  read_csv("./data/analysis/dataset1_oldschool.csv")

# Upload dataset1_staggeredadoption.csv
dataset1_staggeredadoption <-  read_csv("./data/analysis/dataset1_staggeredadoption.csv")

# Upload dataset2_oldschool.csv 
dataset2 <-  read_csv("./data/analysis/dataset2_oldschool.csv") %>% filter(year != 2011)

# Upload dataset2_staggeredadoption.csv
dataset2_staggeredadoption <-  read_csv("./data/analysis/dataset2_staggeredadoption.csv")
########

#### C. Table 2 #### 
m1 <- feols(Tchange_inf2 ~ ND_1.6_bar
            | ADM4_PCODE +year, vcov = ~ADM4_PCODE + year, 
            data = dataset1) 

m2 <- feols(Tchange_inf2 ~ treated
            | ADM4_PCODE +year, vcov = ~ADM4_PCODE + year, 
            data = dataset1_staggeredadoption) 

Bm1 <- feols(npa_after1~ ND_1.6_bar
             | ADM4_PCODE + time, vcov = "twoway",
             data = dataset2)

Bm2 <- feols(npa_after1~ treated
             | ADM4_PCODE + time, vcov = "twoway",
             data = dataset2_staggeredadoption)

etable(m1, m2, Bm1, Bm2, tex = TRUE)
########

#### D. Figure 7 ####
#### 1. Run the regression table first ####
dataset1_staggeredadoption <- dataset1_staggeredadoption %>%
  mutate(nd_mixed = treated*Mixed_inf2,
         nd_encl = treated*Enclave_inf2)

dataset2_staggeredadoption <- dataset2_staggeredadoption %>% 
  mutate(nd_mixed = treated*Mixed_inf2,
         nd_encl = treated*Enclave_inf2)

mod1 <- feols(Tchange_inf2 ~ treated + Mixed_inf2 + Enclave_inf2 +
                nd_mixed + nd_encl 
              | ADM4_PCODE + year, vcov = "twoway",
              data = dataset1_staggeredadoption)

mod2 <- feols(npa_after1~ treated + Mixed_inf2 + Enclave_inf2 +
                nd_mixed + nd_encl
              | ADM4_PCODE + time, vcov = "twoway",
              data = dataset2_staggeredadoption)
########

#### 2. Plot the figure 7a - the shift of territorial control as the DV ####
# Step. 1. Effect of natural disasters among non-mixed and non-enclave setting
mod1_ND <- glht(mod1, linfct = c("treated = 0"))
summary(mod1_ND)
# coef: coef(mod1_ND)
# se: sqrt(diag(vcov(mod1_ND)))
nd <- as.numeric(summary(mod1_ND)$test[[3]]); 

# Step. 2.  Effect of natural disasters among mixed setting 
mod1_ND_mix <- glht(mod1, linfct = c("nd_mixed + treated = 0")); summary(mod1_ND_mix)
## coef: coef(mod1_ND_mix)
## se: sqrt(diag(vcov(mod1_ND_mix)))
nd_mixed <- as.numeric(summary(mod1_ND_mix)$test[[3]]); 

# Step. 3. Effect of natural disasters among enclave setting 
mod1_ND_encl <- glht(mod1, linfct = c("nd_encl + treated = 0")); summary(mod1_ND_encl)
## coef: coef(mod1_ND_encl)
## se: sqrt(diag(vcov(mod1_ND_encl)))
nd_encl <- as.numeric(summary(mod1_ND_encl)$test[[3]]); 


# Step. 4. Create a 5 x 4 matrix to store inputs for creating plot	
plot1 <- as.data.frame(matrix(NA, 3, 3))
colnames(plot1) <- c("estimate", "se", "num")

# Step. 5. Extract the values for creating plot 
plot1[1, 1] <- coef(mod1_ND_encl)
plot1[1, 2] <- sqrt(diag(vcov(mod1_ND_encl)))

plot1[2, 1] <- coef(mod1_ND_mix)
plot1[2, 2] <- sqrt(diag(vcov(mod1_ND_mix)))

plot1[3, 1] <- coef(mod1_ND)
plot1[3, 2] <- sqrt(diag(vcov(mod1_ND)))

plot1$num <- c(3:1)

plot1$estimate <- plot1$estimate * 100
plot1$se <- plot1$se * 100

# Step. 5. Plot it, finally! 
pd <- position_dodge(0.75)

figure7a <- ggplot(plot1, aes(y = estimate, x = num)) + 
  geom_point(position = pd, shape = 19) + 
  geom_linerange(aes(ymin = estimate-1.96*se, ymax = estimate+1.96*se),
                 lwd = 1/2,
                 position = pd) +
  geom_linerange(aes(ymin = estimate-1.645*se, ymax = estimate+1.645*se),
                 lwd = 1,
                 position = pd) +
  geom_hline(yintercept = 0,
             linetype = "dashed") +
  labs(x = "", 
       y = "Change in the likelihood of a shift in territorial control (%)") +
  theme_bw(base_size = 14) +
  annotate("text", y  = 18.91, x = 3.2, label = "18.91") +  
  annotate("text", y  = 3.90, x = 2.2, label = "3.90") +  
  annotate("text", y  = -0.25, x = 1.2, label = "-0.393") + 
  coord_flip() +
  theme(panel.border = element_blank(),
        panel.grid.minor = element_blank(), 
        axis.line = element_line(),
        legend.position = "",
        axis.title=element_text(size=12),
        axis.text.x=element_text(),
        axis.text.y=element_text(),
        aspect.ratio = 0.6)  + 
  scale_x_discrete(limits = c("Homogeneous setting", 
                              "Mixed setting",
                              "Enclave setting"))
########

#### 3. Plot the figure 7b - the battle-related violence as the DV ####
# Step.1. Effect of natural disasters among non-mixed and non-enclave setting 
mod2_ND <- glht(mod2, linfct = c("treated = 0")); summary(mod2_ND)
## coef: coef(mod2_ND)
## se: sqrt(diag(vcov(mod2_ND)))
nd <- as.numeric(summary(mod2_ND)$test[[3]]); 

# Step.1. Effect of natural disasters among mixed setting 
mod2_ND_mix <- glht(mod2, linfct = c("nd_mixed + treated = 0")); summary(mod2_ND_mix)
## coef: coef(mod2_ND_mix)
## se: sqrt(diag(vcov(mod2_ND_mix)))
nd_mixed <- as.numeric(summary(mod2_ND_mix)$test[[3]]); 

# Step.3. Effect of natural disasters among enclave setting 
mod2_ND_encl <- glht(mod2, linfct = c("nd_encl + treated = 0")); summary(mod2_ND_encl)
## coef: coef(mod2_ND_encl)
## se: sqrt(diag(vcov(mod2_ND_encl)))
nd_encl <- as.numeric(summary(mod2_ND_encl)$test[[3]]); 


# Step.4. Create a 5 x 4 matrix to store inputs for creating plot	
plot2 <- as.data.frame(matrix(NA, 3, 3))
colnames(plot2) <- c("estimate", "se", "num")

# Step.5. Extract the values for creating plot 
plot2[1, 1] <- coef(mod2_ND_encl)
plot2[1, 2] <- sqrt(diag(vcov(mod2_ND_encl)))

plot2[2, 1] <- coef(mod2_ND_mix)
plot2[2, 2] <- sqrt(diag(vcov(mod2_ND_mix)))

plot2[3, 1] <- coef(mod2_ND)
plot2[3, 2] <- sqrt(diag(vcov(mod2_ND)))

plot2$num <- c(3:1)

plot2$estimate <- plot2$estimate * 100
plot2$se <- plot2$se * 100


# Step.6. Plot it
figure7b <- ggplot(plot2, aes(y = estimate, x = num)) + 
  geom_point(position = pd, shape = 19) + 
  geom_linerange(aes(ymin = estimate-1.96*se, ymax = estimate+1.96*se),
                 lwd = 1/2,
                 position = pd) +
  geom_linerange(aes(ymin = estimate-1.645*se, ymax = estimate+1.645*se),
                 lwd = 1,
                 position = pd) +
  geom_hline(yintercept = 0,
             linetype = "dashed") +
  labs(x = "", 
       y = "Change in the likelihood of battle-related violence (%)") +
  theme_bw(base_size = 14) +
  annotate("text", y  = 1.14, x = 3.2, label = "1.14") +  
  annotate("text", y  = 3.93, x = 2.2, label = "3.93") +  
  annotate("text", y  = 0.948, x = 1.2, label = "0.948") +
  coord_flip() +
  theme(panel.border = element_blank(),
        panel.grid.minor = element_blank(), 
        axis.line = element_line(),
        legend.position = "",
        axis.title=element_text(size=12),
        axis.text.x=element_text(),
        axis.text.y=element_text(),
        aspect.ratio = 0.6)  + 
  scale_x_discrete(limits = c("Homogeneous setting", 
                              "Mixed setting",
                              "Enclave setting")) 

########

#### 4. Combine the figure 7a and 7b together ####
figure <- grid.arrange(
  figure7a, figure7b)
########
########

#### E. Figure 8 ####
#### 1. Run the regression table first ####
dataset1_staggeredadoption <- dataset1_staggeredadoption %>%
  mutate(nd_mixed = treated*Mixed_inf2,
    nd_encl = treated*Enclave_inf2,
    nd_encl_npa = treated*Enclave_inf2*Influence2,
    nd_mixed_npa = treated*Mixed_inf2*Influence2,
    encl_npa = Enclave_inf2*Influence2,
    mixed_npa = Mixed_inf2*Influence2,
    nd_npa = treated*Influence2)

dataset2_staggeredadoption <- dataset2_staggeredadoption %>%
  mutate(nd_mixed = treated*Mixed_inf2,
    nd_encl = treated*Enclave_inf2,
    nd_encl_npa = treated*Enclave_inf2*Infl2,
    nd_mixed_npa = treated*Mixed_inf2*Infl2,
    encl_npa = Enclave_inf2*Infl2,
    mixed_npa = Mixed_inf2*Infl2,
    nd_npa = treated*Infl2)

mod3 <- feols(Tchange_inf2 ~  treated + Influence2 + nd_npa + Mixed_inf2 + Enclave_inf2 +
                nd_mixed + mixed_npa + nd_npa + nd_encl + encl_npa + nd_mixed_npa + nd_encl_npa
              | ADM4_PCODE + year, vcov = "twoway",
              data = dataset1_staggeredadoption) 


mod4 <- feols(npa_after1 ~  treated + Infl2 + nd_npa + Mixed_inf2 + Enclave_inf2 +
                nd_mixed + mixed_npa + nd_encl + encl_npa  + nd_mixed_npa + nd_encl_npa
              | ADM4_PCODE + time, vcov = "twoway",
              data = dataset2_staggeredadoption) 
########

#### 2. Plot figure 8a - the shift of territorial control as the DV ####
# Step.1. Effect of natural disasters among government-controlled territories in the non-mixed and non-enclave setting 
mod3_homogeneous_gov <- glht(mod3, linfct = c("treated = 0"))
summary(mod3_homogeneous_gov)

# Step.2. Effect of natural disasters among rebel-controlled territories in the non-mixed and non-enclave setting 
mod3_homogeneous_rebel <- glht(mod3, linfct = c("treated + nd_npa = 0"))
summary(mod3_homogeneous_rebel)

# Step.3. Effect of natural disasters among mixed setting controlled by government 
mod3_mixed_gov <- glht(mod3, linfct = c("treated + nd_mixed = 0"))
summary(mod3_mixed_gov)

# Step.4. Effect of natural disasters among mixed setting controlled by NPA 
mod3_mixed_rebel <- glht(mod3, linfct = c("treated + nd_mixed + nd_mixed_npa = 0"))
summary(mod3_mixed_rebel)

# Step.5. Effect of natural disasters among enclave setting controlled by government 
mod3_enclave_gov <- glht(mod3, linfct = c("treated + nd_encl = 0"))
summary(mod3_enclave_gov)

# Step.6. Effect of natural disasters among enclave setting controlled by NPA 
mod3_enclave_rebel <- glht(mod3, linfct = c("treated + nd_encl + nd_encl_npa = 0"))
summary(mod3_enclave_rebel)

# Step.7. Create the dataset with correct coefficient extraction
plot_data <- data.frame(
  Setting = c(
    "Enclave setting (Government)", 
    "Enclave setting (NPA)", 
    "Mixed setting (Government)",
    "Mixed setting (NPA)",
    "Homogeneous setting (Government)",
    "Homogeneous setting (NPA)"),
  estimate = c(
    summary(mod3_enclave_gov)$test$coefficients, 
    summary(mod3_enclave_rebel)$test$coefficients,
    summary(mod3_mixed_gov)$test$coefficients,
    summary(mod3_mixed_rebel)$test$coefficients,
    summary(mod3_homogeneous_gov)$test$coefficients,
    summary(mod3_homogeneous_rebel)$test$coefficients)* 100,
  se = c(
    summary(mod3_enclave_gov)$test$sigma, 
    summary(mod3_enclave_rebel)$test$sigma,
    summary(mod3_mixed_gov)$test$sigma,
    summary(mod3_mixed_rebel)$test$sigma,
    summary(mod3_homogeneous_gov)$test$sigma,
    summary(mod3_homogeneous_rebel)$test$sigma)* 100)

# Step.8. Convert factor levels to ensure Enclave is at the top
plot_data$Setting <- factor(plot_data$Setting, 
                            levels = rev(c("Enclave setting (Government)",
                                           "Enclave setting (NPA)",
                                           "Mixed setting (Government)",
                                           "Mixed setting (NPA)",
                                           "Homogeneous setting (Government)",
                                           "Homogeneous setting (NPA)")))

# Step.9. Define dodge position adjustment for overlapping points
pd <- position_dodge(0.5)

# Step.10. Create the plot
figure8a <- ggplot(plot_data, aes(y = estimate, x = Setting)) + 
  geom_point(position = pd, shape = 19, size = 3) +  # Plot point estimates
  geom_linerange(aes(ymin = estimate - 1.96 * se, ymax = estimate + 1.96 * se),
                 lwd = 1/2, position = pd) +  # 95% CI
  geom_linerange(aes(ymin = estimate - 1.645 * se, ymax = estimate + 1.645 * se),
                 lwd = 1, position = pd) +  # 90% CI
  geom_hline(yintercept = 0, linetype = "dashed") +  # Reference line at 0
  labs(x = "", y = "Change in the likelihood of a shift in a territorial control (%)") + 
  theme_bw(base_size = 14) + 
  # Add text annotations next to each point
  geom_text(aes(label = round(estimate, 2), x = as.numeric(Setting) + 0.3, y = estimate), 
            size = 4, hjust = 0) +  
  coord_flip() +  # Flip axes for better readability
  theme(
    panel.border = element_blank(),
    panel.grid.minor = element_blank(), 
    axis.line = element_line(),
    legend.position = "none",
    axis.title = element_text(size = 12),
    axis.text.x = element_text(),
    axis.text.y = element_text(),
    aspect.ratio = 0.6)

########

#### 3. Plot figure 8b - battle-related violence as the DV ####
# Step.1. Effect of natural disasters among non-mixed and non-enclave setting controlled by the government 
mod4_homogeneous_gov <- glht(mod4, linfct = c("treated = 0"))
summary(mod4_homogeneous_gov)

# Step.2. Effect of natural disasters among non-mixed and non-enclave setting controlled by the rebel 
mod4_homogeneous_rebel <- glht(mod4, linfct = c("treated + nd_npa = 0"))
summary(mod4_homogeneous_rebel)

# Step.3. Effect of natural disasters among mixed setting controlled by government 
mod4_mixed_gov <- glht(mod4, linfct = c("treated + nd_mixed = 0"))
summary(mod4_mixed_gov)

# Step.4. Effect of natural disasters among mixed setting controlled by NPA 
mod4_mixed_rebel <- glht(mod4, linfct = c("treated + nd_mixed + nd_mixed_npa = 0"))
summary(mod4_mixed_rebel)

# Step.5. Effect of natural disasters among enclave setting controlled by government 
mod4_enclave_gov <- glht(mod4, linfct = c("treated + nd_encl = 0"))
summary(mod4_enclave_gov)

# Step.6. Effect of natural disasters among enclave setting controlled by NPA 
mod4_enclave_rebel <- glht(mod4, linfct = c("treated + nd_encl + nd_encl_npa = 0"))
summary(mod4_enclave_rebel)

# Step.7. Create the dataset with correct coefficient extraction
plot_data <- data.frame(
  Setting = c(
    "Enclave setting (Government)", 
    "Enclave setting (NPA)", 
    "Mixed setting (Government)",
    "Mixed setting (NPA)",
    "Homogeneous setting (Government)",
    "Homogeneous setting (NPA)"),
  estimate = c(
    summary(mod4_enclave_gov)$test$coefficients, 
    summary(mod4_enclave_rebel)$test$coefficients,
    summary(mod4_mixed_gov)$test$coefficients,
    summary(mod4_mixed_rebel)$test$coefficients,
    summary(mod4_homogeneous_gov)$test$coefficients,
    summary(mod4_homogeneous_rebel)$test$coefficients)* 100,
  se = c(
    summary(mod4_enclave_gov)$test$sigma, 
    summary(mod4_enclave_rebel)$test$sigma,
    summary(mod4_mixed_gov)$test$sigma,
    summary(mod4_mixed_rebel)$test$sigma,
    summary(mod4_homogeneous_gov)$test$sigma,
    summary(mod4_homogeneous_rebel)$test$sigma)* 100)

# Step.8. Convert factor levels to ensure Enclave is at the top
plot_data$Setting <- factor(plot_data$Setting, 
                            levels = rev(c("Enclave setting (Government)",
                                           "Enclave setting (NPA)",
                                           "Mixed setting (Government)",
                                           "Mixed setting (NPA)",
                                           "Homogeneous setting (Government)",
                                           "Homogeneous setting (NPA)")))

# Step.9. Define dodge position adjustment for overlapping points
pd <- position_dodge(0.5)

# Step.10. Create the plot
figure8b <- ggplot(plot_data, aes(y = estimate, x = Setting)) + 
  geom_point(position = pd, shape = 19, size = 3) +  # Plot point estimates
  geom_linerange(aes(ymin = estimate - 1.96 * se, ymax = estimate + 1.96 * se),
                 lwd = 1/2, position = pd) +  # 95% CI
  geom_linerange(aes(ymin = estimate - 1.645 * se, ymax = estimate + 1.645 * se),
                 lwd = 1, position = pd) +  # 90% CI
  geom_hline(yintercept = 0, linetype = "dashed") +  # Reference line at 0
  labs(x = "", y = "Change in the likelihood of battle-related violence (%)") + 
  theme_bw(base_size = 14) + 
  # Add text annotations next to each point
  geom_text(aes(label = round(estimate, 2), x = as.numeric(Setting) + 0.3, y = estimate), 
            size = 4, hjust = 0) +  
  coord_flip() +  # Flip axes for better readability
  theme(
    panel.border = element_blank(),
    panel.grid.minor = element_blank(), 
    axis.line = element_line(),
    legend.position = "none",
    axis.title = element_text(size = 12),
    axis.text.x = element_text(),
    axis.text.y = element_text(),
    aspect.ratio = 0.6)
########

#### 4. Combine the figure 8a and 8b together ####
figure8 <- grid.arrange(
  figure8a, figure8b)
########
########



